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o 

Space charge effects, being one of the most significant collective effects, play an important role 

(N 

, , in high intensity cyclotrons. Ifowever, for cyclotrons with small turn separation, other existing 

effects are of equal importance. Interactions of radially neighboring bunches are also present, but 
their combined effects has not yet been investigated in any great detail. In this paper, a new 
I I particle in cell based self-consistent numerical simulation model is presented for the first time. The 

model covers neighboring bunch effects and is implemented in the three-dimensional object-oriented 

o 

^ parallel code OPAL-CYCL, a flavor of the OPAL framework. We discuss this model together with 

5/2 its implementation and validation. Simulation results are presented from the PSI 590 MeV Ring 

o 

■ ^ Cyclotron in the context of the ongoing high intensity upgrade program, which aims to provide a 

beam power of 1.8 MW (CW) at the target destination. 

PACS numbers: 29.20.dg;29.27.Bd;41.20.Cv 
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O I. INTRODUCTION 

cn 
o 

Since the invention of the classic cyclotron several decades ago, increasingly higher beam intensities are required to 
^ provide more powerful tools for many new scientific endeavors, such as SpaUation Neutron Sources and in the future 

"x 

;_3 Accelerator Driven Systems that are foreseen to reduce nuclear waste. In high intensity cyclotrons, space charge efi^ects 
play an important role for the following reasons. Firstly, with the absence of longitudinal focusing in cyclotrons, the 
longitudinal space charge force causes additional acceleration for head particles and deceleration for tail particles, 
which may lead to an increase in energy spread. Secondly, there is a strong radial-longitudinal coupling which is 
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influenced by non-linear radial and longitudinal space charge forces. Lastly, the space charge force can reduce the 
vertical tune and increase the vertical beam size, which can cause beam losses when the vertical size extends beyond 
the aperture of the accelerator [1]. 

As shown in [2-7] and experimentally verified at the PSI Injector II [8], an intense particle beam which is properly 
matched to a separate sector cyclotron, develops a spatial stationary circular bunch distribution. Consequently the 
halo is significantly reduced and a corresponding longitudinal decrease of the beam size is observed. We note that 
Kleeven gave in his thesis [9] an "early" hint on this remarkable effect and in [5] a compact derivation can be found. 
On a very practical side, because of the compact stationary distribution we can operate the 3'''^ harmonic resonator, 
the former flat-top resonator in the PSI Injector II, in acceleration mode. 

Non-linear space charge effects in cyclotrons are complex because of the complicated magnetic topology (reference 
trajectory with non constant curvature). Typically there are two approaches to deal with this difficulty: one is an 
approximation with analytic and semi-analytic models which are inexpensive to compute [5] [9] . With such models we 
can obtain a qualitative understanding of the problem. In the past, a number of models have been developed based 
on this philosophy, such as the Sector Model, Disk Model, Sphere Model and the so called Needle Model [3, 4, 10, 11]. 
Another more accurate approach is numerical simulation with macro-particles. In this field, typically two different 
technologies have used to solve space charge fields: Particle-Particle (P-P) methods [12] and Particle-Mesh (P-M) 
methods [6, 7]. In P-P methods, the fields imposed on a given particle are obtained by directly summing up the 
contributions of all other particles at this position. Limited by its low computation efficiency 0{n^) with n denoting 
the number of simulation particles, it is impossible to utilize these methods when the number of particles is large (above 
1 million), even on current state-of-the-art supercomputers. In contrary, in P-M methods the fields are calculated on 
the discrete domains. Due to its high efficiency and high precision, P-M based Particle-In-Cell (PIC) methods [13] are 
widely used in parallel macro-particles simulation codes for different types of accelerators and beam lines [6, 14-17] 
as well as many other areas of computational science, thanks to the development of parallel computation technology 
during recent years. The Parallel PIC model is the method of choice in this study on the beam dynamics of high 
intensity cyclotrons. 

For high intensity cyclotrons, single bunch space charge effects are not the only contribution. Along with the steady 
increase of beam current, the mutual interaction of neighboring bunches in radial direction becomes more and more 
important, especially at large radii where the distances between neighboring bunches diminishes, and even the overlap 
can occur. One example is the PSI 590 MeV Ring Cyclotron [18] with a production current of 2 mA in CW operation 
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and a beam power on target of approximately 1.2 MW. An ambitious upgrade program for the PSI Ring Cyclotron 
is in progress, aimed for 1.8 MW CW beam power on target. The concept involves replacing four aluminum cavities 
by new copper cavities with peak voltages increasing from about 0.7 MV to above 0.9 MV, meanwhile the old flat-top 
cavity remains in use with peak voltage standing at 11.2% of the main voltage. After the planned upgrade is finished, 
the total turn number can be significantly reduced, e.g. from more than 220 turns to less than 170 turns. 

The turn separation is consequently increased as shown in Fig. 1, but remains at the same order of magnitude as the 
measured radial bunch size (at the la level) and is also dependent on the increasing bunch current. Therefore, when 
the beam current increases from 2 mA to 3 mA, the correct treatment of space charge effects is of great importance. 
This includes the mutual space charge effects between radially neighboring bunches. 

Another example is the 100 MeV compact cyclotron (CYCIAE-100) under construction at CIAE [19]. Although its 
beam current is only 200 fiA to 500 /iA, because of the small energy gain per turn (about 200 keV), the turn separation 
is far more smaller than the beam size at outer radii (at the extraction, AR — 1.5 mm) and multiple bunches will 
overlap. 

Because of the complexity of the problem, it is impossible to evaluate neighboring bunch effects precisely and 
self-consistently by explicit analytic expressions. However, high performance computation (HPC) makes it possible 
to treat this problem in greater detail. To our knowledge, little research has been undertaken in neighboring bunch 
effects, and the only published work in that regard was by E. Pozdeyev [7]. He introduced "rigid auxiliary bunches" 
in his serial code CYCO which uses the azimuthal angle as the independent variable. 

In Section II, a new PIC based self-consistent numerical simulation model is presented which, for the first time, 
covers neighboring bunch effects. Section HI describes our three-dimensional object-oriented parallel simulation code 
OPAL-CYCL, a flavor of the OPAL framework. The results of performance testing and code benchmarking are 
presented in Section IV, and followed by its applications on the PSI 590 MeV Ring Cyclotron in Section V. The last 
section is devoted to teh conclusion of the paper. 

II. BASIC EQUATIONS AND PHYSICAL MODEL 
A. PIC MODEL IN CYCLOTRON 

In the cyclotrons and beam lines under consideration, the collision between particles can be neglected because the 
typical bunch densities are low. In time domain, the general equations of motion of charged particle in electromagnetic 
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FIG. 1: (Color) Comparison of calculated turn separation for centroid particles before (red line) and after (blue line) upgrade 
of the PSI Ring Cyclotron. The main cavity voltage is set to 0.735 MV and 0.9 MV for aluminum and copper respectively while 
the flat-top is set to 11.2% of the main voltage. 



fields can be expressed by 



dv{t) 
dt 



q (c;3 X B + E) , 



(1) 



where mo,q,j are rest mass, charge and the relativistic factor. With p — mocyfS we denote the momentum of a 
particle, c is the speed of light, and /3 = {Pxt PyT Pz) is the normalized velocity vector. In general the time [t) and 
position (x) dependent electric and magnetic vector fields are written in abbreviated form as B and E. 
If p is normalized by toqc, Eq. (1) can be written in Cartesian coordinates as 
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The evolution of the beam's distribution function /(x, c/3,t) can be expressed by a coUisionless Vlasov equation: 

df 



(3) 



9t/ + c^-V,/ + g(E + c/3xB).V^^/ = 0, 
where E and B include both external applied fields, space charge fields and other collective effects such as wake fields 



E — Ecxt + Esc, 
B — Bcxt + Bsc- 



(4) 



In order to model a cyclotron, the external electromagnetic fields are given by measurement or by numerical calcula- 
tions. 
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The space charge fields can be obtained by a quasi-static approximation. In this approach, the relative motion of 
the particles is non-relativistic in the beam rest frame, so the self-induced magnetic field is practically absent and the 
electric field can be computed by solving Poisson's equation 

VV(x) = -^, (5) 

where 4> and p are the electrostatic potential and the spatial charge density in the beam rest frame. The electric field 
can then be calculated by 

E = -V(/., (6) 

and back transformed to yield both the electric and the magnetic fields, in the lab frame, required in Eq. (4) by 
means of a Lorentz transformation. Because of the large gap in our cyclotron, the contribution of image charges and 
currents are minor effects compared to space charges [1], and hence it is a good approximation to use open boundary 
conditions. 

The combination of Eq. (3) and Eq. (5) constitutes the Vlasov-Poisson system. In the content followed, the method 
of how to solve these equations in cyclotrons using PIC methods is described in detail. 

Considering that particles propagates spirally outwards in cyclotrons, and the longitudinal orientation changes 
continuously, three right-handed Cartesian coordinate systems are defined, as shown in Fig. 2. The first coordinate 
system is the fixed laboratory frame Slab, in which the external field data is stored and the particles are tracked. Its 
origin is the center of the cyclotron and its X ~ Y plane is the median plane and the positive direction of Z axis 
points to vertical direction. 

The second coordinate system is the local instantaneous frame Siocai, which is a temporal auxiliary frame for the 
space charge solver. Its origin O' is the mass center of the beam and the orientation of the Y' axis is coincident with 
the average longitudinal direction and the positive orientation of the Z' axis points to the vertical direction. 

The third coordinate system is the beam rest frame Sbcam, which is co- moving with the centroid of the beam. It has 
the same orientation and origin as Siocai, but the length in longitudinal direction is scaled by I/7 due to relativistic 
effects. 

At each time step, in order to seek a solution for the space charge fields, the frames Siocai a-nd Sbcam ^^re redefined 
according to current 6D phase space distribution, and all particles are transformed from Slab to Siocai, then a Lorentz 
transformation is performed to transform all particles to Sbcam- The Poisson equation is then solved in the frame 
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FIG. 2: (Color) Schematic plot of the top view of the three coordinates frames. The red curve is the orbit of bunch center, the 
blue area represents bunch shape, and the gray area is the hill region of magnetic field. 



Sbeam- In & 3D Cartesian frame, the solution of the Poisson equation at point {x,y,z) can be expressed by 

1 



G{x, x' , y, y' , z, z')p{x\ y' , z')dx' dy dz , 



(7) 



with G the 3D Green function 



G{x,x',y,y',z,z') 



(8) 



^{x-x')^ + {y-y')^ + {z-z') 
assuming open boundary conditions. The typical steps of calculating space charge fields using the Hockney's FFT 
algorithm [13] is sketched in Algorithm 1, where the quantities with superscript D (discrete) refer to grid quantities. 

With respect to the external magnetic field two possible situations can be considered: in the first situation, the 
real field map is available on the median plane of the existing cyclotron machine using measurement equipment. In 
most cases concerning cyclotrons, the vertical field, B^, is measured on the median plane (z — 0) only. Since the 
magnetic field outside the median plane is required to compute trajectories with z ^ 0, the field needs to be expanded 
in the Z direction. According to the approach given by Gordon and Taivassalo [20], by using a magnetic potential 
and measured on the median plane at the point (r, 6, z) in cylindrical polar coordinates, the 3rd order field can 
be written as 

OB, 1 



Br{r,e,z) = z^~^Z^Cr, 

or D 

Bg{r,e,z) = ~ p—'^e, 

r ot) or 

B,{r,e,z) = B,- iz^C 



(9) 



7 



Algorithm 1 3D Space Charge Calculation 

1: procedure 3DSpaceCliarge(In: p, G, Out: E^cBsc) 

2: Create 3D rectangular grid wliich contains all particles 

3: Interpolate the charge q of each macro-particle to nearby mesh points to obtain 

4: Lorentz transformation to obtain in the beam rest frame Sboam 

5: FFT p° and G° to obtain p° and G° 

6: Determine on the grid using = • G^ 

7: Use FFT^i of 0° to obtain (jy^ 

8: Compute E° = -V<;/>^ 

9: Interpolate E at the particle positions x from E^ 

10: Inverse Lorentz transform to obtain Esc and Bsc in frame Siocai and transform back to Siab 

11: end procedure 
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All the partial differential coefficients are computed on the median plane data by interpolation, using Lagrange's 
5-point formula. 

In the other situation, 3D field for the region of interest is calculated numerically by building a 3D model using 
commercial software during the design phase of a new cyclotron. In this case the calculated field will be more accurate, 
especially at large distances from the median plane i.e. a full 3D field map can be calculated. For all calculations in 
this paper, we use the method by Gordon and Taivassalo [20]. 

Finally both the external fields and space charge fields are used to track particles for one time step using a Ath 
order Runge-Kutta (RK) integrator, in which the fields are evaluated for four times in each time step. Space charge 
fields are assumed to be constant during one time step, because their variation is typically much slower than that of 
external fields. 
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B. NEIGHBORING BUNCH EFFECTS 

The code is intended to model steady state conditions for the muhi-bunch beam dynamics. In cyclotrons the pattern 
of turn separation AR is affected by many factors. These include machine characteristics such as the magnetic field, 
the acceleration voltage profile, the accelerating phase of the RF resonators and initial centering conditions of the 
injected bunches. Generally, in cyclotrons, AR reduces gradually with increasing beam energy. For machines like the 
PSI Injector II, AR stays sufficiently large from injection to extraction, and in such cases, neighboring bunch effects 
are negligible. For others, like the PSI 590 MeV Ring Cyclotron under consideration in this section, AR decreases 
strongly during the course of acceleration, which results in the need of considering the neighboring bunch effects in 
order to obtain a precise description of the beam dynamics. In our model, we apply an iterative procedure to determine 
the number of bunches necessary for a converged simulation. Initially a single bunch with phase space density /q and 
the average radial position R^ is injected. This single bunch is tracked with space charge for one revolution period 
Tr- Then the new average radial position of the bunch R^ and the bunch rms size Trms — \/ ^rms Vms ^'^^ calculated 
from the actual particle distribution. The turn separation AR at injection position is then given by AR = R^ — Rs- 
If the condition: 

AR<M X rrms, (11) 

is fulfilled (where M is a parameter given by the user) , the 6D phase space is stored as fn^ . The code is switched to 
multi-bunch mode, and fji^ will be used as the initial phase space for the following (Nb — 1) neighbouring bunches 
which will be injected one by one per Tr time, where Nb is the number of neighboring bunches given by the user. 

If the condition of Eq.(ll) is not fulfilled, the value of R^ is assigned to the variable Rg. This single bunch is tracked 
with space charge for another revolution period T^, and the new average radial position of the bunch Re, the bunch 
rms size Tims and the turn separation AR are calculated again. After that, the condition in Eq.(ll) is evaluated and 
the same procedure is repeated accordingly. 

The underlying assumption for this Ansatz is that all bunches have the same phase space distribution when they 
reach a certain position, i.e. when fji^ is saved. This is realistic and reasonable when the machine is running in 
a steady state. It need be mentioned that up to that position the coherent instabilities which might be caused by 
neighboring bunch effects, are not covered by this model. 

This procedure is summarized in Algorithm 2. 

In the multi-bunch algorithm above, two parameters M and Nb are introduced to set the time of injecting new 
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Algorithm 2 Multi-Bunch Injection Algorithm 



1: procedure Injection(In: TVs, M, fo, Rs) 

2: Inject /() at radius Rs, total bunch number ijvs ~ 1 

3: Track fo for one revolution period Tr, obtain fa^ 

4: Calculate radius Re and bunch size rrms 

5: while {Re — Rs > M x rrms) do 

6: Save Re — s> Rs 

7: Track (r^ for Tr, re-calculate Re and rrms 

8: end while 

9: Save fn^ 

10: while (ijvs < Nb) do 

11: Inject /r^ and increment ijvs 

12: Track all bunches for Tr 

13: end while 

14: end procedure 




FIG. 3: (Color) Schematic plot of the top view of 5 bunches and the grid of computation domain. The grid size on X' — Y' 
plane is NxxNy, and the broken lines represent the orbits of bunch centers. 



bunches and the total bunch number respectively. The proper settings of these two parameters are crucial for the 
precise evaluation of neighboring bunches effects. 

In order to quantify the range of the two parameter Nb and Ad, let's consider a 2D non-relativistic DC beam. The 
Bassetti-Erskine [21] formula for the electric field of a 2D Gaussian charge distribution is in general an analytical 
expression in terms of the complex error function. In case of an axisymmetric and Gaussian charge distribution the 
electric field can be expressed by 



Esc(r) 



27reo/3 cr 



Ur, r < a 



r > a. 



(12) 
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FIG. 4: (Color) The radial space charge field Esc{r) of the multiple Gaussian beams at the central beam region. The beams 
are truncated at 10% of their charge density and rms beam size is r^ms = 1.0 (mm). All graphs are on the same scale. 

In this expression, a is the truncated radius, Iq the beam current and unit vector in radial direction. Using Eq.(12), 
it is easy to calculate the electrostatic field generated by Nb Gaussian beams (with Nb an odd number) which are 
all on a straight line. We can now study the effect of Nb and M on the centre beam Nbc — {Nb + l)/2 of the 
configuration as shown for several configurations in Fig. 4. 

Although in the cyclotron, the situation is much more complex (the charge distribution and radial position of the 
bunches are time dependent quantities), we use this model to obtain approximate initial conditions for iV^ and M 
and hence can estimate the contributions for bunches away from Nbc- 

For instance, the setting with Nb — 9 and M = 4.5 gives convergent results for the PSI Ring Cyclotron with 3 mA 
beam current. 

In theory, when the maximum bunch number Nb equals to the total turn number of the machine, one can eventually 
obtain the fully self-consistent solution of the problem within our model. In reality, it is not feasible to simulate a full 
set of bunches, which typically range from several ten to several hundred. The scale of the number of particles and 
the dimensions of the needed grid are still beyond the capability of today's supercomputer resources. 

In a multi-bunch simulation the energy of the bunches at different turns can be substantially different. For a 
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multi-bunch simulation with 9 neighbouring bunches, if the kinetic energy of the first bunch is Ek = 100 MeV and 
energy gain per turn is AE^ = 2 MeV, hence the velocity difference of the first and last bunch is 6.5%. In this case 
there is no single rest frame in which the relative motions of particles are non-relativistic, as required by our scheme 
to calculate the space charge forces. Consequently it is not sufficient to use only one rest frame and one single Lorentz 
transformation. In order to calculate the space charge fields more precisely, we use an adaptive binning technique 
outlined in Ref. [22] (section IV). We note that in our application neither radiation nor retardation effects play a 
significant role and can therefore be neglected. In the rest frame of the beam, transverse currents effects can also be 
neglected and hence no longitudinal magnetic field component must be considered. 

First we create the same amount of energy bins (Nb) as we have bunches in the simulation. An average relativistic 
factor 7i for the i'^ bunch with iV* simulation particles is computed. 



Then every particle is grouped into the energy bin whose is closest to its 7. In this way, the energy spread of each 
bin is small and the relative motions of the particles in the same bin are small. After binning we perform the Lorentz 
transformation, calculate the space charge field and perform back-transformation for each bin respectively. Finally 
the field data is summed up to give the total space charge force imposed on each particle. 

The energy spread of the bunches can be large (MeV range), especially in cyclotrons without flat-top cavities, and 
at large radii. This may result in an overlap of energy distributions of neighbouring bunches, and hence the energy 
bins have to be recalculated i.e. all particles need to be regrouped after each revolution. It is worth noting that, in 
cyclotrons, the energy difference of neighboring bunches changes with the increasing radius. Therefore the energy 
difference of the neighboring bins is not constant. Specifically, the energy difference between the i*^ and (i— 1)**^ bins, 
AEi^i^i, differs with the energy difference between the {i + 1)"^ and the i^^ bins AiJ^+i^i. 



The above model and algorithm are implemented in the object-oriented parallel PIC code OPAL-CYCL. OPAL- 
CYCL is one of the flavors of the OPAL (Object Oriented Parallel Accelerator Library) framework [23]. This framework 
is a powerful tool for charged-particle optics in general accelerator structures and beam lines using the MAD languages 
with extensions. OPAL is based on the CLASSIC [24] hbrary and the IP^L framework [25]. The CLASSIC library 
is a C-I--I- class library which provides services for building portable accelerator models and algorithms and inputing 
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III. IMPLEMENTATION WITHIN THE OPAL FRAMEWORK 
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language to specify complicated accelerator systems in general. IP^L is an Object-Oriented C-|— I- class library which 
provides abstractions for particles and structured field calculation in a data-parallel programming style. It provides 
an integrated, layered system of objects. The upper layers contain global data objects of physical/mathematical 
quantities, such as particles, fields and matrices of meshes and typical methods performed on these objects such as 
differential operators and multi dimensional FFT's. The lower layers contain the objects relevant to parallelization 
such as data distribution, domain decomposition, communication among processors, load balancing and expression 
templates. Statistical data, such as root mean square (rms) quantities, are computed on the fly (in situ) and stored 
in conjunction with phase space and field data in the H5Part [26] file format. In a post processing step the data can 
be analyzed using the visualization tool HSPartROOT [27]. 

In addition, apart from the multi-particle simulation mode, OPAL-CYCL also has two other serial tracking modes 
for conventional cyclotron machine design. One mode is the single particle tracking mode, which is a useful tool for 
the preliminary design of a new cyclotron. It allows to compute basic parameters, such as reference orbit, phase shift 
history, stable region and matching phase ellipse. The other one is the tune calculation mode, which can be used to 
compute the betatron oscillation frequency and Vz- This is useful for evaluating the focusing characteristics of a 
given magnetic field map. 

A more detailed description of the hierarchical layout, the parallelization and the implementation issues of the 
OPAL framework and OPAL-CYCL code can be found in the User's Reference Guide [23]. 

IV. PERFORMANCE TEST AND VALIDATION 

In order to evaluate the performance and to benchmark the functionalities of the newly developed code, we performed 
different types of simulations on the 72 MeV Injector II cyclotron of PSI, which has been intensively studied before. 
Some results are presented in this section. 

A. Single particle tracking and tune calculation 

In theory, there is an eigen-ellipse for any given energy in a cyclotron under stable conditions. When the initial phase 
space matches this eigen-ellipse, the oscillation of the beam envelope amplitude will be minimal and the transmission 
efficiency will be maximal. In the present test, the eigen-ellipse at 2 MeV kinetic energy is calculated using the single 
particle tracking mode of OPAL-CYCL. The result is compared to FIXPO [28, 29]. At PSI the FIXPO code has been 
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FIG. 5: (Color) Radial eigen-ellipse of 2MeV at the symmetric line of a magnet in Injector II cyclotron. 
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FIG. 6: (Color) Tune diagram of Injector II cyclotron, compared with FIXPO code. 

the standard simulation tool for the design and parameter optimisation of the Injector II and the Ring Cyclotron as 
well as for a gas driven muon trap in a cyclotron shaped magnetic field. The code integrates single particle orbits by 
use of a predictor corrector algorithm up to the third order. In Fig. 5 the matched radial ellipse with an initial offset 
of X = 2.0 mm, px = O.Omrad at the symmetry line of the sector field is shown. Excellent agreement is obtained when 
the time step is set to 1 ps in OPAL-CYCL, although FIXPO uses a different tracking algorithm with the azimuthal 
angle as the independent variable. 

The tune diagram of Injector II is computed using the tune calculation mode of OPAL-CYCL, as shown in Fig. 6. 
The result from FIXPO and OPAL-CYCL are again in good agreement even though different numerical algorithms 
are used. 

The field interpolation scheme, particle tracking and tuning calculation functionalities are validated substantially 
by the above tests. 




14 



10 





1 


1 1 

Total time - 


— 1 






Self fields ■■ 


- -X 






Integration •■ 





- 




Data dump 


B - 


55 






)(-..,"■■"• 










. )K. 








"■■■X.,., 


,., % 








""^ SK-... 
















X-. 








B 


l 






B 









16 32 64 128 



Number of processors 



FIG. 7: (Color) The time consumption as a function of processors on Cray XT3, CSCS. 



B. Parallel scalability test 



In order to observe the parallel scalability of the code, we have performed a detailed study of strong scaling, i.e. 
the problem size remains constant while increasing the number of computing resources. One million particles are 
used and tracked 200 time steps on the Injector II. The initial beam has a Gaussian type distribution. The grid size 
is 64 X 64 X 64 which is decomposed onto a two dimensional grid of processors. All the intermediate phase space 
data is dumped into a single H5Part file. The dynamic load balancing frequency as well as the phase space dumping 
frequency are set to 10. The results are shown in Fig. 7. 

We can see that good scalability is achieved up to 128 processors. Above 128 processors, the time consumption 
of the phase space dumping starts to become significant. The reason for the behavior is the increasing overhead in 
communication with respect to the amount of data to be stored Nevertheless, the scalability of the space charge solver 
and the particles integrator still benefit from a large number of processors. 



C. Stationary Round Distribution in the PSI Injector II 

Space charge effects usually result in an increase in beam size and emittance. That is detrimental to beam dynamics. 
However there are cases, where space charge effects can actually play a positive role. The PSI Injector II cyclotron is 
a space charge dominated machine, in which a very compact stable beam is developed within the first several turns, 
and thereafter, the charge distribution does not change significantly. This stationary situation remains essentially 
unchanged until extraction and the beam phase width is about 2° in the last turn. This is due to the combined effect 
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of the strong coupling between the radial and longitudinal directions in the cyclotron and the space charge when 
the beam current increases above 1mA. S. Koscielniak and S.Adam reproduced this phenomenon by using the serial 
two-dimensional code PICN [4] . PICN is based on the Needle model, which treats the beam as an ensemble of charged 
vertical needles with the same height as the beam. In this model, the vertical motion of particles is separated from 
the horizontal component and the internal motion within needles is neglected. In order to validate the space charge 
solver module of our code, we performed the 1 mA, 3 MeV coasting beam simulation on this machine and compared 
the resuh with PICN. 

We used the same initial distribution as in Ref. 30. Sciongitudinai = 13.4mm ( 15° phase width), 2(Ttransvcrsc = 
2.52 mm. The initial emittances of the radial and azimuthal directions are set to zero. In the vertical direction, 
2(72 = 4 mm and 2(7^^ = 3.68 mrad. The total macro-particles number is 1 million. Figure 8 shows the top view of 
beam shape in the local beam frame. We can see that a stable core is developed after about 10 turns, which is faster 
than the formation of stable haloes. 

When comparing Fig. 8 with Fig. 3 and Fig. 4 in Ref. 30 calculated by PICN, we can find the results agree with 
each other qualitatively. Both of these two codes predict the formation of a compact stable core and wide haloes after 
40 turns. However, there still exist visible differences. PICN shows an almost "round" charge density distribution, in 
the case of OPAL-CYCL we still see low density halo. The longitudinal size of the core is about 2 mm longer than 
the transverse size. This difference mainly comes from different physical models used in the two codes. PICN uses a 
so called smooth approximation which treats the particle orbits as pure circles, but with a sinusoidal radial/vertical 
focusing(betatran oscillation) having a realistic value of about 1.14 at 3 MeV. Of course, no realistic magnetic field 
can be azimuthally constant and in the same time focusing in both directions [31]. In addition all particles in close 
proximity to the horizontal plane are represented by a single "needle"; in OPAL-CYCL only those particles close to 
each other in configuration space are represented by a single macro-particle. The latter is believed to be more realistic 
and accurate. 

It can be concluded from this comparison that OPAL-CYCL can reliably reproduce the stable "round" beam 
formation caused by space-charge effects in Injector II and hence can accurately reproduce the single bunch dynamics 
in cyclotrons. 
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FIG. 8: (Color) Top view of a 1mA, 3MeV coasting beam in PSI Injector II in the local frame Siocai of PSI Injector II. Up: 
turn 0, 5, 10. Down: turn 20, 30, 40. To compare with figures in Ref. 30, the beam's transport direction is along the negative 
direction of the abscissa axis. 

V. APPLICATIONS 

We start this section by describing the two cyclotrons under consideration. The key parameters of the machines 
are given in Table I. In addition we mention the fundamental RF frequency is 50.633 MHz. In the Injector II due to 
the circular bunch formation, discussed in section IV, the original designed third harmonic flat-top resonator is now 
being used as an additional accelerating structure. It is now obvious to ask the question if one can find a feasible 
working point for the Ring Cyclotron with the same characteristics as obtained in the Injector II. However because 
of the overlapping tuns in the Ring Cyclotron, the situation is much more complex than in the Injector II. Those two 
issues are addressed for the first time in the remainder of this paper. 
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TABLE I: Key parameters of the two sector cyclotrons 





orbit radius kin. energy 


avg. power 


cLvg. field 


peak field ma,^ 


^n. rigidity harmonic number resonators / 




(m) (MeV) 


(MW) 


(T) 


(T) 


(Tm) cavities 


Injector II 


0.44. ..3.3 0.87. ..72 


0.15 


0.33. ..0.36 


1.08 


1.25 10 4 


Ring Cyclotron 


2.1... 4.45 72. ..592 


1.3 


0.6. ..0.9 


2.17 


4.0 6 4 (1 flat-top) 



A. Different phase width studies of the PSI Ring Cyclotron 



Although a very compact beam with a phase width of about 2° can be extracted from the Injector II, it is nevertheless 
subject to the expansion in the longitudinal direction in the 72 MeV beam transfer line because of space charge effects 
and chromatic dispersion. For the future 3 mA beam, this will have a significant impact on the beam dynamic of 
the Ring Cyclotron. In response to this, a rebuncher running on the 10th harmonic is planned to be installed on the 
beam line to make bunches as short as possible at the injection point of the Ring Cyclotron. The final bunch length 
achieved is a critical aspect of the Ring Cyclotron. 

In order to obtain a clear perspective on this issue, OPAL-CYCL was applied to do numerical simulation by tracking 
Gaussian type beams with 3 different initial conditions. The initial longitudinal phase widths (Gtr) are set to 2°, 6° 
and 10°, respectively, and the initial energy spread is neglected. The initial conditions of the horizontal and vertical 
directions are identical. The initial distribution is not correlated in phase space. The simulation used 10^ macro- 
particles and 32x32x32 gird sizes. The peak voltage of the four main resonators and the 3rd harmonic flat-top 
resonator are 0.9 MV and 0.403 MV (11.2% of accelerating voltage), respectively. The time step is set to 0.1 ns. It 
takes about 7 hours on CRAY XT3 of CSCS using 64 processors to track particles from the injection to the extraction. 

Figure 9 shows the development of the beam rms size on the transverse and the longitudinal direction. We can 
see the beam is compressed gradually in the longitudinal direction. Meanwhile, in the transverse direction, the beam 
size increases fast during the first several turns because of the mismatch of initial conditions. Thereafter the beam 
size does not change significantly until the beam arrives at the extraction region where it is distorted by the external 
magnetic field (extraction bump). Figure 10 shows the projection of phase space onto the mid plane of the machine, 
and Fig. 11 plots the histogram along the longitudinal direction at 112° azimuthal position of turn 0, 50 and 150. We 
can see for the bunch with the initial phase width of 2°, the bunch maintains a very compact shape with a stable 
round core without haloes. When the initial phase width increases, the size of the core only widens slightly (less than 
5 mm), while the spiral tails expand in the longitudinal direction and are unable to develop stable haloes. However, 
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FIG. 9: (Color) Comparison of the rms beam size in the transverse direction (left) and longitudinal direction (right) at 112° 
azimuthal position of each turn in PSI Ring Cyclotron. 



the beam does not expand notably in the radial direction, which means no substantial increase of the beam loss on 
the extraction septum is expected for bunches with initial phase width less than 10°. 



B. Neighboring bunch effects in the PSI Ring 

As discussed in Section I, neighboring bunch effects may have an appreciable influence on beam dynamics in the 
Ring. This can be evaluated by comparing the difference in single bunch and multiple bunch simulations as described 
in Section III. We have run simulations for 3, 5, 7 and 9 bunches and found that the difference between the 7 bunch 
scenario and that of 9 bunch scenarios is small, i.e. is viewed as converged, as illustrated in Fig. 12 and discussed 
already in Section II B. From Fig. 13 we conclude that the FWHM of the transverse profile is reduced by approximately 
33% comparing a 1 and 9 bunch simulation. For the energy spread (FWHM) we have a reduction in the order of 14% 
i.e from 0.7 MeV to 0.6 MeV. 

As Gordon explained in [10] the particle motion in a cyclotron is always perpendicular to the force resulting in a 
vortex motion. In order to obtain the observed sharpening of the distribution we need an additional, azimuthal force. 
A possible explanation of the origin of this force is due to the observed broken circular symmetry when considering 
neighboring bunches in the simulation. However more efforts are needed in order to understand this effect in greater 
detail. 

From the comparison we conclude that the integration of neighboring bunch effects into the model has non-negligible 
impacts on the beam dynamics for beam currents beyond 1 mA in the PSI Ring. The bunch becomes more compact 
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FIG. 10: (Color) Top view of 3 mA bunch distributions with 2°, 6° and 10° initial phase widths at the initial position(top) 
turn 50 (middle), and 150 (bottom) in the local frame Siocai of 112° azimuthal position of PSI Ring Cyclotron. 
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FIG. 11: (Color) Histogram of 3mA bunch distributions with 2°, 6° and 10° initial phase widths at initial posit ion(left), turn 
50 (middle), and 150 (right) in PSI Ring Cyclotron. 
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FIG. 12: (Color) Top view of 1mA bunch distributions at the turn 130 in the local frame Siocai at 112° azimuthal position of 
turn 130 in PSI Ring Cyclotron. The results are obtained from single bunch (left), 7 bunches (middle) and 9 bunches (right) 
simulations, respectively. 
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FIG. 13: (Color) Comparison of the histograms along the transversal direction in the local frame Siocai (left) and the energy 
spectra (right) of 1mA beam at 112° azimuthal position of turn 130 in PSI Ring Cyclotron. 

in the transverse direction and the energy spread is slightly reduced. Therefore neighboring bunch effects have a 
positive influence on reducing beam loss in high intensity operation. 



VI. CONCLUSIONS AND DISCUSSIONS 



A physical model for the beam dynamics in high intensity cyclotrons, which includes for the first time the space 
charge effects of neighboring bunches, is presented in this paper. This model is implemented in an object-oriented 
three-dimensional parallel PIC code (OPAL-CYCl), as a flavor of the OPAL framework. 

The performance tests on the CRAY XT3, CSCS demonstrate a good scalability of OPAL-CYCL with respect to 
the number of used processors. The three operation modes of this code (tune calculation, single and multiple particle 
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mode) are validated by code comparison. OPAL-cycl has been successfully applied to study the behavior of the PSI 
Ring Cyclotron at high intensities. 

The beam intensity of this high power facihty is practically hmited by uncontrolled losses at the extraction element 
of the cyclotron, originating from beam tails that are developed during the acceleration process. As shown in the 
results, the generation of beam tails can be avoided if short bunches with a phase length of 2° or less are injected. An 
upgrade plan is under way to generate such short bunches with the help of a 10th harmonic buncher. Furthermore 
it is observed that the neighboring bunch effects can help to narrow the transverse beam size and reduce the energy 
spread. The differences between single and multi-bunch simulations are in the order of 33% and 15% in the compared 
quantities, beam size and energy spread respectively. This is a significant difference between single and multi-bunch 
simulation, and hence justifies the presented simulation method and their application. This is an important step 
towards the quantitative understanding of beam tails in high power cyclotrons. Considering the fact that in the PSI 
Ring Cyclotron the total sum of controlled and uncontrolled losses are in the order of 10~^, a precise beam dynamics 
simulation must cover radial neighboring tuns in order to predict losses with the mentioned intensities. 

It is planned to refine these simulations within the next few years by a more detailed determination of the initial 
particle distribution at the injection of the PSI Ring Cyclotron. A quantitative comparison of the results with 
measured beam quantities will be presented in a future paper. 
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